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1 Introduction 

Discontinuous Galerkin (DG) methods offer an enormous flexibility regarding local 
grid refinement and variation of polynomial degrees rendering such concepts pow- 
erful discretization tools which have proven to be well-suited for a variety of differ- 
ent problem classes. While initially the main focus has been on transport problems 
like hyperbolic conservation laws, interest has meanwhile shifted towards diffusion 
problems. Specifically, we focus here on the efficient solution of the linear systems 
of equations that arise from the Symmetric Interior Penalty DG method appUed 
to elUptic boundary value problems. [IJ The principal objective is to develop ro- 
bust preconditioners for the full "DG-flexibility" which means to obtain uniformly 
bounded condition numbers for locally refined meshes and arbitrarily (subject to 
mild grading conditions) varying polynomial degrees at the expense of Unearly scal- 
ing computational work. A first step towards that goal has been made in [3] treating 
the case of geometrically conforming meshes but arbitrarily large variable polyno- 
mial degrees which already exposes major principal obstructions. In this paper we 
complement this work by detailed studies of several issues arising in [3]. 

To our knowledge the only concept yielding full robustness with respect to poly- 
nomial degrees is based on Legendre-Gaufi-Lobatto (LGL) quadrature. Specifically, 
in the framework of auxiliary space methods low order finite element discretiza- 
tions on LGL-grids can be used to precondition high order polynomial discretiza- 
tions. However, when dealing with variable degrees the possible non-matching of 
such grids at element interfaces turns out to severely obstruct in general the de- 
sign of efficient preconditioners. To overcome these difficulties we propose in [3] a 
concatenation of auxiUary space preconditioners. In the first stage the spectral DG 
formulation (SE-DG) is transferred to a spectral continuous Galerkin formulation 
(SE-CG). In the second stage we proceed from here to a finite element formulation 
on a specific dyadic grid (DFE-CG) which is associated with an LGL-grid but be- 
longs to a nested hierarchy. The latter problem can then be tackled by a multilevel 
wavelet preconditioner presented in forthcoming work. The overall path of our it- 
erated auxiliary space preconditioner therfore is SE-DG SE-CG DFE-CG. It 
should be noted that a natural alternative is to combine the first stage with a domain 
decomposition substructuring preconditioner as proposed in [6] admitting a mild 
growth of condition numbers with respect to the polynomial degree. 
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We are content here for most part of the paper with brief pointers to the detailed 
analysis in [3], [4] and [2] to an extent needed for the present discussion. 

Section 2 introduces our model problem, the LGL technique is explained in Sec- 
tion 3. The auxiUary space method is detailed in Section 4, while Sections 5 and 6 
consider stages 1 and 2 of our preconditioner. Finally in Section 7 we give some 
numerical experiments that shed hght on the constants that arise in four basic in- 
equalities used in the second stage. 

2 Model problem and Discontinuous Galerkin formulation 

Given a bounded Lipschitz domain £2 CW' with piecewise smooth boundary we 
consider as a simple model problem the weak formulation: find u £ Hq{Q) such 
that 



of Poisson's equation —Au = / on X2 with zero Dirichlet boundary conditions m = 
on dQ. For simplicity, we assume that Q is the union of a collection ^ of finitely 
many (hyper-)rectangles, which at most overlap with their boundaries. More com- 
plex geometries can be handled by isoparametric mappings. By .^i{R) we denote 
the /-dimensional facets of a (hyper- )rectangle R and by = Ur^^^i {R) the union 
of all these objects. Let Hk (R) be the side length of R in the k-th coordinate direction. 

The polynomial degrees used in each cell R are defined as p = {pk)i=it where 
is the polynomial degree in the k-\h coordinate direction. We introduce the piecewise 
constant function 5 = {H,p) that collects the hp approximation parameters. On 5 
we impose mild grading conditions, see [3] for the details. 

For the spectral discretization of our model problem, we choose the DG spectral 
ansatz space Vg '.^ {v e L^{Q) : v\r e QpiR) for all R e where Qp{R) is the 
tensor space of all polynomials of degree at most p on the (hyper-)rectangle R. 

We employ the standard notation of DG methods for jumps and averages on the 
mesh skeleton and on dQ . 

The Symmetric Interior Penalty Discontinuous Galerkin method (SIPG) ag (u,v) = 
(/, v) for all V € Vg with the SIPG bilinear form 



as{u5,V5):= ^ (VM5,Vv5)i;+ ^ (-({VM5},[v5])F-([M5],{Vv5})f) 

+ L 7(Of{\u5\ , [vsDf = {f,V5)Q, vs e Vs 




with (Of := max {(Of R- , (Of R+ } for internal faces F and (Of R± : 
For boundary faces FcdQwe set (Of^R := 
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3 Legendre-GauB-Lobatto (LGL) grids 

Denoting by {^i)f^i the zeros of the first derivative of the p-th Legendre polyno- 
mial Lp, (in ascending order), and setting = —I and = I, '^p = {^i)o<i<p 
is the Legendre-GauB-Lobatto (LGL) grid of degree p on the reference interval 
/ = [— 1 , 1], see e.g. [5]. In combination with appropriate LGL weights (w,)o<Kp the 
LGL points of order p can be interpreted as quadrature points of a quadrature rule 
of exactness 2p—l. In [4] we prove quasi-uniformity of the LGL-grids {^p)peNf 
i.e., remains uniformly bounded independent of p, where hj = — ^i-i \ for 
l<i<p-l. 

The particular relevance of tensor product LGL-grids for preconditioners for 
spectral element discretizations lies in the two norm equivalences (see [5]) 

\MHi{R)^\\^M\H'(R) for^l <peQp(/?), iG{0,l}, (1) 

which hold uniformly for any J-dimensional hypercube R= xf^ilk where J^^^ is 
the piecewise multi-Unear interpolant on the tensor product LGL-grid. 

4 Abstract theory: Auxiliary Space Method 

The auxiliary space method (ASM) [9, 11, 10] is a powerful concept for the con- 
struction of preconditioners that can be derived from the fictitious space lemma 
[8,7, 9]. 

Given a problem a{u,v) — f{v) for all v G V on the linear space V equipped 
with a bilinear form a(-, •) : V x V ^ M, we seek an auxiliary space V with an 
auxiliary form •) : V x V ^ M that is in some sense close to the original one but 
easier to solve. Note that we neither require V C V nor V C V which is important 
in the context of non-conforming discretizations. Therefore on the sum V — V + V 
we need in general another version fl(-,-):yxy— >IRas well as a second form 
b{-,-) :V xV ^M. which dominates a on V and plays the role of a smoother. The 
required closeness of the spaces V and V is described with the aid of two linear 
operators Q :V -^V and Q .V ^V. Specifically, these operators have to satisfy 
certain direct estimates involving the above biUnear forms. For the details on the 
ASM conditions see [9]. 

Lemma 1 (Stable Splitting [9]). Under the assumptions of the ASM, we have the 
following stable splitting 

a(v,v)~ inf {b{w,w) +d{v,v)) forallvGV. 

weV,veV: v=w+Qv 

The main result of the ASM is given in the following theorem [9]. 

Theorem 1 (Auxiliary Space Method). Let Cb and be symmetric precondi- 
tioners for B and A, respectively. Let S be the representation ofQ:V^V. Then 
Ca := Cb -I- SC^^^S-^ is a symmetric preconditioner for A. Moreover, there exists a 
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uniform constant C such that the spectral condition number ofC\A satisfies 
ff(CAA)<Cff(CBB)jf(C^A). 

For a given practical application it remains to identify a suitable auxiliary space 
V, the bilinear forms d:V xV and d,b :V xV ^R, as well as the two linear 
operators Q and Q, such that ASM conditions are satisfied. In addition efficient 
preconditioners for the "easier" auxiUary problems and Cb need to be devised. 
Of course, the rationale is that the complexity to apply and Cr should be much 
lower than solving the original problem. 

Note that the operator Q need not be implemented but enters only the analysis. 

5 Stage 1: ASM DG-SEM ^ CG-SEM 

In the first stage, we choose the largest conforming subspace V : = n //q (i2 ) of 
V := Vs as auxiliary space so that Q can be taken as the canonical injection. The 
definition of the operator Q can be found in [3]. 

The main issue in this stage is the choice of the auxiliary form b{-,-). Using LGL- 
quadrature combined with an inverse estimate for the partial derivatives in the form 
a{-,-) we arrive at 



b{u,v):= £ £ u{^)v{^)c^W^, W^:= 




Here the weights 1 are chosen as 

^ ' \ J3ici, else, 

where Wp ji± is the LGL quadrature weight on F seen as a face of and the param- 
eters j3,pi can be used to "tune" the scheme. The resulting mattix B is diagonal so 
that the application of Cb := requires only £>'{N) operations. It is shown in [3] 
that all ASM conditions are satisfied for this choice of ft ( ■ , ■ ) . Numerical experiments 
show that the parameters j3i and pi can by and large be optimized independently of 
the polynomial degrees. 

6 Stage 2: CG-SEM ^ CG-DFEM 

The second stage involves three major ingredients, namely 

(1) the choice of spaces of piecewise multi-linear finite elements on hierarchies of 
nested anisotropic dyadic grids, to permit a subsequent application of efficient 
multilevel preconditioners, 

(2) the construction of the operators Q and Q, and 

(3) the construction of the auxiliary bilinear form b{-,-). 
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As for (1), the non-matching of LGL-grids for different degrees p at interfaces pre- 
vents us from taking low order finite element spaces as auxiliary space for the high 
order conforming problem resulting from the first stage. Instead, with each LGL- 
grid '^p we associate a dyadic grid ^D,p, which is roughly generated as follows: 
starting with the boundary points { — 1,1} as initial guess we adaptively refine the 
grid. A subinterval in the grid is bisected into two parts of equal size, if the small- 
est of the overlapping LGL-subintervals is longer than a times its length. The pa- 
rameter a therefore controls the mesh size of the dyadic grid. However, for input 
LGL-grids of different polynomial degrees the resulting dyadic grids are not nec- 
essarily nested yet. How to ensure nestedness while keeping the grid size under 
control is shown in [3]. The key quality of the associated dyadic grids ^o.p is that 
mutual low order piecewise multi-hnear interpolation between the low order finite 
element spaces on 'S^p{R) ,^D,p{R) is uniformly //^-stable, see [3] for the proofs. 
Denoting by Vk,D,p{R) the space of piecewise multi-linear conforming finite ele- 
ments on 'i^D,p{R), we now take V := Vg nH^{Q) and V := Vhj)nH^{n), where 
Vh,D = {v e C^iQ) ■.yRG^,viR := VR G Vh.oA'^)}- 

Concerning (2), the operator Q is defined element-wise as follows. For a given 
element vertex z S =^o(^) let p* denote the polynomial degree vector whose kih 
entry is the minimum of the Mi entries of all degree vectors associated with elements 

sharing z as a vertex. Here a grading of the degrees is important. Let <P^ G Qi (R) 
the multi-hnear shape function on R satisfying conditions ^z(y) = dy^^ for all y £ 
^o{R)- Then, we define 

K-=^h,D,Pti^zVR)&VH,D,pt{R) and v* = >|v*eQ^jW, (2) 

where ■J^hDp*'^p* are the dyadic piecewise multilinear and high order LGL- 
interpolants on the respective grids. Summing-up over the vertices of R, we define 

Vr:= L n^yh,D,p{R) and Qrvr := v*r := £ v,* € Qp(/?) . (3) 

ze.^oiR) ze^o(R) 

The operator Q is defined analogously with the roles of dyadic and LGL-grids ex- 
changed, see [3]. 

To finally address (3), for the structure of the form /?(•, •) from the first stage the 
direct estimates in the ASM conditions are no longer valid. It has to be suitably 
relaxed along the following lines. We make an ansatz of the form 

Kv,w):=i:i:( b%^S,{^,w)+ £ 4,.a(v,w)), (4) 

where 3^o,k{R) is the collection of those LGL-subcells Sc, £ E xf^^{l,. . .,pk{R)} 
with side lengths hf'^ in the LGL-grid that are strongly anisotropic according 

to (max/^^/z|^'')//i[^*^ > Caspect for a fixed constant Caspect > 0> while S^i^k{R) is 
comprised of the remaining "isotropic" cells. On the isotropic cells in =^ we 
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use an inverse estimate applied to piecewise multi-linear LGL-interpolants of v and 
w. On the remaining anisotropic cells we retain integrals over the variable involving 
the partial derivative and use quadrature in the remaining variables. For this auxiUary 
form b{-,-) and the above operators Q and Q we can verify all ASM conditions, see 
[3]. Note that the Gramian B is no longer diagonal and we refer to [3] for efficient 
reaUzations of Cb. 

7 Numerical experiments: Constants in the basic interpolation 
inequalities 

A fundamental role in the proof of the ASM-conditions in the second stage SE-CG 
DFE-CG is played by four basic interpolation estimates. In particular, knowing the 
size of the constants arising in these inequalities and their dependence on the poly- 
nomial degrees helps understanding the quantitative effects observed in more com- 
plex situations later on. 

As before, let 0^ denote the affine shape function now on the reference interval 
/ = [— 1, 1] C M satisfying <P^{x) = 5j forx,z G {—1, 1}. By .i^^, we denote the poly- 
nomial interpolation operator on the LGL-grid "i^g for polynomial degree q and by 
^h,D,q the piecewise affine interpolation operator on the dyadic grid ^D,q associated 
with 

A major tool for proving the ASM conditions is given by the following theorem. 
Theorem 2. Assume that cp <q< p for some fixed constant c > 0. Then we have 

\'^q{^zV)\H'niT)^MH'n{I) /o?" «// V € Qp(/), Z € {-1, 1}, W € {0, 1}, (5) 

and 

\^h,D,qi^zv)\H'n{r)SMH'"it) /o'-a«v€y/,,D,p(/), ze{-i,i}, me{o,i}. (6) 

We determine next numerically the smallest constants that fulfill the inequalities 
(5) and (6). This can be obtained by solving generalized eigenvalue problems for the 
largest generalized eigenvalue. For all dyadic grids we choose the grid generation 
parameter a = 1.2, which balances two effects: on the one hand, the generated 
auxiliary space is rich enough for a good approximation while on the other hand, to 
keep the solution of the auxiliary space feasible, the dyadic grid does not have too 
many degrees of freedom. Figiu-e 1 shows the dependence of the smallest possible 
constants on the polynomial degrees p and q in the range I < p,q <64. 

We observe that the constants in (5) and (6) become large for m = when the 
quotient p/q increases, but eventually stay bounded as long as cp < q < p for a 
fixed c > 0. For m = 1 we find uniform moderate constants in (5) and (6) for ar- 
bitrary choices of p and q. While the nodes in the LGL-grids move gradually with 
increasing degree the associated dyadic grids change more abruptly which explains 
the jumps in the graph in Figure 1(c). 

We are particularly interested in the behavior of the constants when the quo- 
tient of p and q is fixed, i.e., we restrict ourselves to a cross section through the 
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(c) (6), m = (d) (6), m = 1 

Fig. 1 Dependance of the constants in (5) and (6) on p and q. 

3-dimensional plots along a line in the p^-plane. As an example, we choose p = 2q 
representing strongly varying degrees on adjacent elements. The smallest constants 
in the inequalities for polynomial degrees q up to 128 are displayed in Figure 2. 




Fig. 2 Constants in the basic interpolation inequalities for p = 2q (dashed line: (5), solid line: (6)). 



While for m — the constants quickly approach an asymptotic value for both 
(5) and for (6), this is not true for (5) and m = 1. In this case we observe a very 
slow monotonic convergence to its asymptotic limit. Thus for moderate polynomial 
degrees one still observes a significant growth. Since this estimate is relevant for the 
ASM conditions on the operator Q in the second stage, this leads to some growth 
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of the condition number of the preconditioned problem for moderate polynomial 
degrees and significant inter-element jumps, although it eventually stays uniformly 
bounded independent of the polynomial degree q. 

8 Summary and outlook 

In this paper we sketch a preconditioner for the spectral symmetric interior penalty 
discontinuous Galerkin method that, under mild grading conditions, is robust in 
variably arbitrarily large polynomial degrees, announcing detailed results given in 
[3]. The concept is based on the LGL-techniques for spectral methods combined 
with judiciously chosen nested dyadic grids through an iterated application of the 
auxiliary space method. A detailed exposition of a multiwavelet preconditioner for 
the dyadic grid problem, an extension to locally refined grids with hanging nodes, 
strategies for parallel implementations, and the treatment of jumping coefficients 
will be presented in forthcoming work. 
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